% This script generates ensembles of temperature trajectories for figure 3. The same ensembles are used for table 1 in the SI, and for the estimation of the change in the SCC in the SI. 
close all
clear all
 
% Load RCP data (This script is run once for each RCP. Each time the script is run, you must "comment out" all but one of the RCPs below.)
%load rf_rcp8_5.txt
%load rf_rcp6_0.txt
%load rf_rcp4_5.txt
load rf_rcp2_6.txt
%yrs=rf_rcp8_5(:,1);
%yrs=rf_rcp6_0(:,1);
%yrs=rf_rcp4_5(:,1);
yrs=rf_rcp2_6(:,1);
%forcing=rf_rcp8_5(:,2);
%forcing=rf_rcp6_0(:,2);
%forcing=rf_rcp4_5(:,2);
forcing=rf_rcp2_6(:,2);

% Set start time in 2020
tstart=256;

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Specify initial temperature anomaly in 2020
DeltaT0 =1; % degrees C

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Specify model parameters %%%%%%%%%%%%%%%%%%
% We will measure time in years,
% Specify size of time interval between points.
% so need a factor of number of seconds in year
% Effectively divides time in seconds by yearlength
dt=1; 
% So then have to multiply speed by yearlength
yearlength=3.1 *1e7;

% Specify number of points at which solution is wanted
nPeriods=length(forcing(tstart:end)); % Same as RCP dataset

% Specify # of intermediate steps taken between the sample points
nSteps=12; % Monthly

% Number of temperature trajectories without noise
Ntrials0=1;

% Number of temperature trajectories with noise
Ntrials=1000;

% Effective heat capacity
C=1*1e9; % Watts/metre^2 /Kelvin 

% Standard deviation of noise
sigmaQ=.9375e8;

% Range of values of equilibrium climate sensitivity and feedback parameter
ECS = 0.5:0.25:20
lambdas = 3.7 ./ ECS

parfor k = 1:1:length(lambdas)
    CS = ECS(k);
    Sigma=0;
    lambda=lambdas(k);   % Joules/metres^2/Kelvin 
    Speed = lambda*yearlength/C; % Kelvin per year 

    % And convert F/lambda to a function
    Level=ts2func(forcing(tstart:end)/lambda);

    [DeltaTnoiseless,Times0]=Hasselmann(Speed,Level,Sigma,DeltaT0,Ntrials0,nPeriods,nSteps,dt);

    Sigma=sigmaQ/C;
    [DeltaT,Times]=Hasselmann(Speed,Level,Sigma,DeltaT0,Ntrials,nPeriods,nSteps,dt);

	%Create a CSV write file
	% First suppress the unit 3rd dimension in Delta T matrix
	Temp=squeeze(DeltaT);
	% Then concatenate
	M=[Times DeltaTnoiseless Temp];
 
	% filename = ['RCP85output_' num2str(CS) '.txt']; % Write output to these files for RCP8.5
	% filename = ['RCP60output_' num2str(CS) '.txt']; % Write to these files for RCP6.0
	% filename = ['RCP45output_' num2str(CS) '.txt']; % Write to these files for RCP4.5
	filename = ['RCP26output_' num2str(CS) '.txt']; % Write to these files for RCP2.6
	csvwrite(filename,M)

end